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Abstract 


A semi-Lagrangian method for parabolic problems is proposed, 
that extends previous work by the authors to achieve a fully conser¬ 
vative, flux-form discretization of linear and nonlinear diffusion equa¬ 
tions. A basic consistency and convergence analysis are proposed. 
Numerical examples validate the proposed method and display its 
potential for consistent semi-Lagrangian discretization of advection- 
diffusion and nonlinear parabolic problems. 
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1 Introduction. 


It is well known that, in their original formulation, Semi-Lagrangian 
(SL) schemes were developed for the case of purely hyperbolic equa¬ 
tions and do not guarantee conservation of mass. In recent years, an 
increasing number of SL discretization approaches have been proposed 
to circumvent these limitations. 

Various extensions of the SL strategy to second order, parabolic 
problems, appear for example in a, m, esi, m, n m (see ai so a 
for a more complete review), as opposite to the “classical” technique 
of treating in Eulerian form diffusion terms within SL schemes. In 
spite of the variety of applications involved, the common feature of 
these works is to replace the computation of the solution at the foot 
of a characteristic with a weighted average of values computed at mul¬ 
tiple points. A rigorous derivation of these techniques can be based 
on the Feynman-Kac representation formula, with stochastic trajec¬ 
tories integrated backward in time as in conventional SL schemes. 
However, due to their stochastic origin, all these generalizations of 
the SL approach are based on the trace form of parabolic problems, 
and unsuitable for the divergence form 

u t = div (DVu) (1) 

which is the most widely used in many applications, especially those 
to computational fluid dynamics on geophysical scales, for which SL 
methods for advective problems have proven to be especially useful 
(see e.g. the review in [3]). In the previous paper [2], we have proposed 
for the first time a SL method for parabolic problems in divergence 
form. Numerical experiments reported therein show that the method, 
albeit formally only first-order accurate in time, allows to compute 
remarkably accurate approximations of linear and nonlinear parabolic 
problems, as well as to achieve easily higher order accuracy in space 
also on nonuniform meshes. 

On the other hand, the last two decades have witnessed an increas¬ 
ing interest in conservative, flux-form SL (FFSL) schemes which have 
been proposed, e.g., in [21], na, Ha, ng, m, eh, esj, eh, ebi, for 
various linear and nonlinear advective problems. These methods are 
exactly mass conservative, as opposed to methods based on upstream 
remapping of computational cells (also called inherently conservative 
in the literature), see e.g. [5], [6], (12[. 11]. [16], [IT. [22] . [26]. 

EH, [32], m, M- 

The method proposed in [2], however, even though discretizing 
Equation ([Tj) directly, is not formulated in flux form at the discrete 
level and does not guarantee exact mass conservation. In the present 
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work, we try to fill the gap by extending the approach of [2] to achieve 
a fully conservative, flux form discretization of the diffusion equa¬ 
tion ([Tj) . Along the same lines as [2], we outline a consistency and 
stability analysis of the method in a simplified setting, and perform a 
number of numerical simulations to validate the proposed method. 

The paper is structured as follows. In Section [2j our novel FFSL 
discretization is introduced, while an analysis of its consistency and 
stability properties will be presented in Sections [3] and [4j Section [5] 
treats the extensions of the scheme to problems in higher dimensions. 
Finally, numerical results obtained with the proposed method for both 
linear and nonlinear models will be presented in Section [6] while some 
conclusions on the potential advantages of our approach will be drawn 
in Section 0 


2 Flux Form Semi-Lagrangian meth¬ 
ods for second order problems. 

In order to sketch the basic ideas of the scheme, we restrict for the 
moment to the approximation of the pure diffusion equation 

ut = (v(x,t)u x ) x (2) 

in a single space dimension, considered for simplicity on the whole real 
line. The extension to advection-diffusion equations and to multiple 
dimensions will be discussed later in this Section and in Section [5l 
We have shown in [2] that © can be approximated via an abstract 
difference operator in the form of a convex combination of pointwise 
values. More specifically, let At and Ax denote the time and space dis¬ 
cretization steps, respectively, with t n = nAt for n E [0,T/At]. The 
space grid is supposed to be infinite and uniform, that is, for i £ Z, 
Xi = iAx. In the case of the conservative scheme, we will also con¬ 
sider the intermediate points xy- 1/2 = (i ± 1/2) Ax, which will appear 
as endpoints of grid cells 12* = [xj_i/ 2 , ^+ 1 / 2 ] • The nonconservative 
method is derived by setting 

u(xi,t n+ 1 ) « Afu(xi + 6f,tn) + A~u(xi - 8~,t n ), (3) 

provided the constants Af, A~, 8f. 8f satisfy the consistency condi¬ 
tions 

' A t + A 7 = i + o(At 2 ) 

Af8+-A-5- = Atu x + 0(At 2 ) 

* A+8f + Ar6f = 2Atu + 0(At 2 ) 

At8f-AT5f = 0{At>). 
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Note that, actually, the additional condition df 4 = 0 (At 2 ) should also 
be added, which is however already implied by the second condition. 
We will denote the numerical solutions at time t n by the vector V n = 
( u f)«ez- While in the nonconservative form SL scheme these values 
should be understood as pointwise values, in a FFSL scheme they 
must be interpreted as cell averages. 

In [2], the abstract structure introduced in Q takes the specific 
form of the nonconservative approximation 

<*' = lynfe + >t )+lynfe - v>. w 

where the pointwise values are reconstructed by the interpolation op¬ 
erator (of degree p ) I p [V n ], while the displacements Sf are given as 
solution of the equation 

Sf = \J2At v{xi ± 5 f,t n ). (6) 

For a flux form variant of the above approach, rather than approxi¬ 
mating a pointwise value, the cell averages 

Uj_(t) « - 1 - [ u(x,t)dx 
Ax J n . 

of the solutions of Equation © must be approximated, as customary 
in finite volume methods, to obtain the discrete mass balance 

Ui(t n + 1) ~ Ui(t n ) + 1/2 — 1/2) • CD 

The flux might be defined in abstract form as 

\ / f-xk+Sk r%k \ 

= o ( / u(x,t n )dx — / u(x, t n )dx ) 

* \Jx k Jx k -5 k J 

for a properly defined 4- The practical version of the scheme imple¬ 
ments this operator with a numerically reconstructed function, so that 
the proposed FFSL method for Equation Q is defined as 

< +1 =<+h (*r + i /2 - (s) 

with a numerical flux given, for k = i A 1/2, by 

1 / fx-k+Sk rxk \ 

F k=~( R q [V n ](x)dx - / R q [V n )(x)dx , (9) 

z \Jx k Jx k —5 k J 

and with 

4 = \j 2 At v(x k ,t n ). ( 10 ) 
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Here, the operator i2 9 [V] is a polynomial reconstruction of degree q 
which satisfies over each cell Q m = [x m _i/ 2 , ® m +i/ 2 ] the properties 
of having the correct cell average w m and reconstructing with order 
0{Ax q+1 ) a smooth solution. There is a close relationship between the 
two operators I p and R q . We refer to pH] for an in-depth discussion of 
these theoretical issues, while in Section [4] will recall the main points of 
interest for the purpose of a stability analysis of the proposed method. 
It should be remarked that, in this case, it is not required to solve any 
equation to obtain 5 k . Furthermore, notice that both methods ([5]) 
and ([8]) employ viscosity values frozen at the time level t n , so that an 
extension to the nonlinear case is straightforward. 

In order to provide a more intuitive interpretation of the proposed 
scheme, it could be observed that the effect of diffusion (as modeled 
by Fick’s law implicit in Equation ([2])) is to move mass from regions 
with high density towards regions with low density. The difference of 
integrals in ([9]) can be shown to be an approximate model of this pro¬ 
cess. Indeed, the interval y/2At v is the correct space scale associated 
to the given diffusion over a time step At. Furthermore, it can be 
observed that 

[ x k+5k / \ 

J u(x, t n )dx ~ 6 k u + —,t n J , 



6 k u 



so that, as a consequence, one obtains 


™ „ h ( ( 5 k \ ( 5 k \\ 

?k ~ -j l« l*fc + -j ,tn ) - u { Xk ~ ~J ,tn )) 

5 2 

~ Y U x (x k ,t n ) = A tv(x k ,t n ) u x (x k ,t n ) 

which clearly makes Q consistent with ([2]). Of course, this is rather 
a heuristic explanation than a rigorous consistency analysis. Consis¬ 
tency (and its optimal rate) will be proven in the next Section on the 
basis of conditions ©■ 

We now show briefly how the proposed FFSL method can be com¬ 
bined with analogous methods for the flux form advection equation. 
Consider the advection-diffusion equation 


u t + ( f(x,t)u) x = (v{x,t)u x ) x , 


( 11 ) 


taken again for simplicity on the infinite real line. The combination 
of the FFSL methods for advection and diffusion is obtained, for the 
purposes of this paper, by a simple operator splitting approach. This 
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combination results in a method that is first order accurate in time, 
which is compatible with the time accuracy of the FFSL method for 
diffusion, as it will be seen in Section [3j In a first step, any FFSL 
method for advection can be used. In Section [6| we consider for con¬ 
creteness the well-known method presented in p~5i . This results in 
a numerical approximation of the cell averages of the solution 
of with v = 0. An approximation of the solution of the complete 
Equation ( [XT] ) is then obtained by computation of formula ([8]) based 
on the intermediate values 1 . 


3 Consistency. 

To prove consistency of the flux form SL method, we start by neglect¬ 
ing the error associated to the space discretization and rewriting the 
right-hand side of ([7]) in the more convenient form: 


u i(tn ) + (-^i+l /2 ^h-1/2) — u i(t n ) + 


2Ax 


c i + l/2 + 'h + l/2 


C i+1/2 


r i+l/2 


' x i+l/2~$i+l/2 


r x i-l/2+$i-l/2 f x i- 1/2 

J + J 

Jx i- 1/2 Jx 


— Ui (fn) T 


u(x , t n )dx = 

x i+l/2 — <h+l/2 f x i+1/2 


i—1/2 —1/2 , 


+ 


L 


2Ax \Jx i _ 1/2 -s i _ 1/2 


"i+l/2+<h+l/2 


i-l/2+<h-l/2 


rx 

•> x i- 1/2 

r x i + l/2 \ 

/ w(x,t n )dx. 

d x i — l/2 J 


+ 


We obtain therefore 
1 


Ui{tn)+ Ax (^ +1 /2 ^- 1 / 2 ) - 2A.t ( J x 


x i+l/2~ <h+l/2 r x i+l/2+3i+l/2 \ 

+ / u(x,t n )dx. 

Z ^ X ' ' x i-l/2~^i-l/2 J x i-l/2+Si-l/2 J 

( 12 ) 


The terms dj-i-i /2 can also be expressed as 
4tl/2 = \J 2 At v{Xi±i/2) = 


(13) 


1 2At ( ^ ± + ^~v X x + O (Ax 3 ) ) = 

1 2At v ( 1 ± _ + — + o (Ax 3 ) = 


2u 8 v 


. Ax v x Ax 

( 1 i ^ ^ l VxT. 


4 1 / 


I61/ 


2zy 2 


+ O (Ax 3 )^ , 
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where, in the last row, we have applied the Taylor expansion of the 
square root and collected all terms of order 0( Arc 3 ). It is now easy to 
recognize that, once replaced (up to a 0( Ax 2 ) error) the integrals in 
the last row of (12) with their midpoint approximation, i.e., 

i+l/2±&i+l/2 


L 


u(x,t n )dx = (Ax±5 i+1 / 2 : f5 i _ 1 /2)u(xi±di) + 0 (Ax 2 ) 


x i-l/2^i-l/2 

and using (13), then (12) is in the form ([3]), with 
St = 5i 

l 6 

Ax ± (6j + 1/2 — <5j_ 1 / 2 ) 


4 * = 


2Ax 


(14) 

(15) 


= 


We check now that (14)- (16) satisfy the consistency conditions Q. 
The first condition is clearly satisfied. As for the second, we have 

A+St-A-S- = S,(A*-A-) = 

= A tu x + O (Ax 2 ) . 

The third condition is rewritten as 

Afdf + A-Sf = 5 2 (At + A~) = 

= 2At v (l + O (Ax 2 )) . 

Finally, for the fourth condition we obtain 

A+8f-Arsf = (2A^) 3 / 2 ^ + 0(Ax 2 ) 

= O (At 2 ) , 


in which the last expression is achieved by taking into account that 

Therefore, the abstract operator is consistent. Finally, introducing 
the polynomial reconstruction R q [V], we note that fluxes are approx¬ 
imated with accuracy of order 0(Ax q+2 ). Then, the estimate 

L(Ax,At) = o(V + ^Q, (16) 


holds for the truncation error L, with r = 1 and s = q + 2. 













4 Stability. 

We briefly discuss the stability of the proposed scheme in the case of 
a constant coefficient equations, while variable coefficient equations 
seem to require a deeper study. First, we note that a well understood 
case is the construction of the nonconservative scheme with I p in the 
form of a symmetric Lagrange interpolation, for p odd. For exam¬ 
ple, if p = 3 , the value R\V]{x) is computed, for x G [x k ,Xk+i], by 
interpolating the values v k -i, ■ ■ ■ ,v k + 2 - In the pure advection case, 
this construction is known to be stable (see 0) In order to con¬ 
struct a stable conservative scheme, we can define the reconstruction 
R q (see m for more details and examples) according to the axioms 

i) For x G fl m , R q [W](x) = Q m (x) for some polynomial Q m G P ? , 
with q even; 

ii) For x G fl m , Q m {%) depends on the values w k , with k = m — 
q/2,... ,m + q/2, and more precisely it satisfies the conditions 

} / Q m (x)dx = w k ( k = m- q/2, + q/2). 

Jn k 


With this definition, it is proven in |10] that, if p = q + 1, then 
^ rx-\- Ax/2 


Ax 


' x—Ax/2 


R q [W}(OdZ = I P [W] (x) 


(17) 


Consider now the problem (12) for a constant diffusivity u. In this 
case 6i = 5, so that, rewriting (121) for the reconstructed numerical 
solution and using ©, we get 


v? +1 = v? + 


1 




-FT 


~Ax \ i+l > 2 ~ '- 1 / 2 


(18) 


l r*i+i/ 2 -S i r 

JX;_t/o—5 Jx 


x i-l/2 


E i+l/2+<5\ 

r i- 1 / 2+5 J 


R q [V n ]{x)dx = 


= - I p [V n ](x-S) + - I P [V n ](x + 6). 

The scheme can then be recast as the convex combination 

V n+1 = -B + V n + -B~ V n , 

2 2 


(19) 


where the terms B±V n represent respectively the schemes written in 
extended form as 

v? +1 = I p [V n ]( Xi ±6), 
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which (see [8]) are stable in the 2-norm for any 5 and p. This implies 
that 

ll B± ll 2 < 1, 

and therefore that the complete scheme is also stable in the same 
norm. 


5 Multiple space dimensions. 

In this Section, we discuss the extension of the proposed approach to 
the d-dimensional case. The extension is straightforward in the case 
of a structured orthogonal grid and diagonal diffusivity matrix 


A (x,t) = diag(zq (x,t ),..., u d (x, t)), (20) 

and in particular, for a variable but isotropic diffusion (for which 
u\(x,t) = ■ ■ ■ = Vd(x,t)). Then, the diffusion equation reads 


ut = div(A(x, t)Vu) = 

d 

3 = 1 

By a first-order expansion of the solution with respect to time, we get 


u(x,t + At) = 


i(x,t) + At^2(vjU Xj ) Xj (x,t) + 0(At 2 ) = 

3 = 1 


E 

3 = 1 


u(x, t) At, 1 . . 

+ — (ai'jUx, ) Xj (x, t) 


d d 
1 d 


+ 0(At 2 ) = 


= u[ 


(x,t) + 2^2 [A t(dvjU Xj ) Xj (x,t)} + 0(At 2 


j =i 


which shows that, up to hrst-order accuracy, the d- dimensional version 
can be obtained by averaging the diffusion operators in each direction, 
with the only modification of scaling each one-dimensional diffusivity 
by a factor d. This allows to split the diffusion along each of the 
variables, where (in case of an orthogonal mesh) flux tubes are also 
aligned with the grid. 

For example, the 2-dimensional equation 


u t = (zq (x,t)u Xl ) xl + (. v 2 (x,t)u X2 ) X2 
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would be approximated by the scheme (written with an obvious nota¬ 
tion at the node Xi = (a^,®^)) 


v? +1 =v? + 


- ft 


/± x 2 ^11+1/2,12 r h-l/2,i 2 + -^1,12+1/2 ^it,12+1/2J > 

( 21 ) 

in which, for example, the flux F ™_, 2 i2 would be given by 


- FJ' 


pn 

r ii+l/2,i 2 


1 

4 


E i 2 + l/2 / -a: i 1 + l/2 + ^i 1 +l/2,i2 


c i 2 -l/2 J3: i 1 +l/2 


r x i 2 +i/2 n 
J Xi 


, X i „_ 1 /2 Jx il+1 /2~ $i 1 +l/2, i 2 


R[V n ](x)dx- (22) 
i2[F n ](x)dx ) , 


and the displacement 5j 1+1 /2,2 2 defined, at the centers of E and W 
facets of a square cell, by 


<5*1+1/2,22 = ^/4At i/i(x il+1/2 ,i 2 ,i n ). 


Note that the integrals appearing in all the fluxes of the form (22) 
have an integration domain consisting of a row (or column) of grid 
cells, among which at most one is intersected. 


6 Numerical experiments. 

Several numerical experiments have been carried out with simple im¬ 
plementations of the FFSL method proposed above, in order to assess 
its accuracy and stability features also in cases more complex than 
those allowing a complete theoretical analysis. The accuracy of the 
proposed discretization has been evaluated against analytic solutions 
or reference solutions obtained by alternative discretizations in space 
and time. Due to the close relationship between the methods proposed 
here and those in [2], the choice of the test cases follows closely the 
outline in our previous paper, in order to allow for a clear comparison 
between the conservative and nonconservative SL discretizations. 

6.1 Constant coefficient case. 

In a first set of numerical experiments, the constant coefficient diffu¬ 
sion equation 

u t = vu xx x G [0, L\ 

was considered, on an interval [0, L] with L = 10. Periodic boundary 
conditions were assumed and a Gaussian profile centered at L /2 was 
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Resolution 

/ 2 rel. error (SL) 

Z 2 rel. error (FFSL) 

N 

M 

h 

h 

h 

Ro 

r 2 

200 

50 

0.8 

1.50 • 10" 2 

3.33 • 10~ 4 

1.5 • 10" 2 

3.79 • 10~ 4 

200 

100 

0.4 

1.45 • 10~ 2 

2 . 9 i • io ~ 4 

1.45 • 10~ 2 

2.02 • 10 -4 

200 

200 

0.2 

6.53 • 10" 2 

5.89 • 10 -4 

6.53 • nr 2 

2.28 • 10" 4 

400 

100 

1.6 

6.53 • 10~ 3 

i .92 • nr 4 

6.53 • nr 3 

1.78 • 10 -4 

400 

200 

0.8 

1.48 • 10~ 2 

8.29 • nr 5 

1.48 • 10~ 2 

9.45 • 10” 5 

400 

400 

0.4 

1.43 • 10" 2 

7.43 • 10~ 5 

1.43 • 10~ 2 

5.03 • nr 5 


Tabic 1: Relative errors for the constant coefficient pure diffusion case in the 
Z 2 -norm, nonconservative (SL) and conservative (FFSL) scheme, first and 
third order space discretizations. 


considered as the initial condition. In this case, the exact solution 
can be computed up to machine accuracy by separation of variables 
and computation of its Fourier coefficients on a discrete mesh of N 
points with spacing Ax = L/N. We consider the FFSL method de¬ 
scribed in the previous Sections on a time interval [0, T] with T = 2, 
with time steps defined as At = T/M. The stability parameter of 
standard explicit discretizations of the diffusion operator is defined as 
/i = uAt/Ax 2 . We consider the case with v = 0.05 first, whose rela¬ 
tive errors are reported in Tables 00 as computed m the / 2 and Iqq , 
respectively. The parallel results for the advection diffusion case 

u t + au x = vu xx x G [0, L\ 

with a = 1.5, v = 0.05 are reported in Tables [3][4j respectively. In this 
case, the Courant number is defined as C = a At/Ax and the advective 
flux was computed using either linear or quadratic reconstruction, 
along the lines of the well known method m- The fractional flux 
employed yields in the low Courant number case a discretization that 
is equivalent to a flux form Lax-Wendroff method. 

It can be observed that, in the linear case, the errors obtained by 
the SL and FFSL methods are in general of the same order of mag¬ 
nitude, with consistently smaller values for the conservative variant. 
The slightly larger errors of the FFSL case in the advection diffusion 
cases are partly explained by the fact that, in the present implementa¬ 
tion, for simplicity a second order reconstruction was employed for the 
advective terms, compared to the third order accurate reconstruction 
employed for the diffusion terms. 
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Resolution 

Zoo rel. error (SL) 

Zoo rel. error (FFSL) 

N 

M 

h 

h 

h 

Rq 

r 2 

200 

50 

0.8 

1.73 • 10" 2 

3.75 • 10~ 4 

1.73 • 10” 2 

4 ^ 

L- 

4 ^ 

M 

O 

1 

it. 

200 

100 

0.4 

1.67-10~ 2 

6.46 • 10~ 4 

1.67- 10" 2 

2.36 • 10~ 4 

200 

200 

0.2 

7.41 • 10" 2 

1.63 • HT 3 

7.41 • 10" 2 

2.67- 10“ 4 

400 

100 

1.6 

7.56 • 10~ 3 

2.85 • 10~ 4 

7.56 • 10~ 3 

2.08 • 10~ 4 

400 

200 

0.8 

1.71 • 10“ 2 

8.86 • 10” 5 

1.70 • 10" 2 

1.11 • 10” 4 

400 

400 

0.4 

1.65 • 10" 2 

1.94 • 10~ 4 

Or 

O 

1 

to 

5.89 • 10~ 5 


Table 2: Relative errors for the constant coefficient pure diffusion case in the 
Zoo-norm, nonconservative (SL) and conservative (FFSL) scheme, first and 
third order space discretizations. 


Resolution 

Z 2 rel. error (SL) 

Z 2 rel. error (FFSL) 

N 

M 

C 

h 

h 

h 

Rq 

R‘2 

200 

50 

1.2 

0.8 

1.21 • 10” 2 

3.77-10- 4 

2.67 • 10 -2 

6.46 • 10~ 4 

200 

100 

0.6 

0.4 

3.38- 10“ 2 

2.91 • 10~ 4 

4.83 • 10" 2 

1.89 • 10” 3 

200 

200 

0.3 

0.2 

4.15 • 10“ 2 

1.57-10~ 3 

h— 1 

h-L 

O 

1 

2.67- 10“ 3 

400 

100 

1.2 

1.6 

5.02 • 10” 3 

1.84 • 10” 4 

1.25 • 10“ 2 

2.14 • 10” 4 

400 

200 

0.6 

0.8 

1.30 • 10~ 2 

1.16 • 10 -4 

3.21 • 10~ 2 

4.79 • 10" 4 

400 

400 

0.3 

0.4 

2.95 • 10" 2 

4.27-10~ 4 

4.41 • 10" 2 

6.67- 10" 4 


Table 3: Relative errors for the constant coefficient advection diffusion case 
in the / 2 -norm, nonconservative (SL) and conservative (FFSL) scheme, first 
and third order space discretizations. 


Resolution 

^OO 

rel. error (SL) 

Zoo rel- error (FFSL) 

N 

M 

c 

h 

I 

1 

h 

Ro 

R-2 

200 

50 

1.2 

0.8 

1.39- 

10" 2 

4.85 • 10~ 4 

3.05 • 10" 2 

6.06 • 10“ 4 

200 

100 

0.6 

0.4 

3.87- 

10" 2 

2.91 • 10 -4 

5.50 • 10" 2 

1.91 • 10~ 3 

200 

200 

0.3 

0.2 

4.73- 

10" 2 

1.71 • 10~ 3 

1.28 • 10" 1 

2.69 • 10~ 3 

400 

100 

1.2 

1.6 

5.82- 

10~ 3 

2.56 • 10" 4 

1.45 • 10~ 2 

2.39 • 10" 4 

400 

200 

0.6 

0.8 

1.49 • 

10“ 2 

2.51 • 10~ 4 

3.68 • 10- 2 

4.85 • 10~ 4 

400 

400 

0.3 

0.4 

3.38- 

10“ 2 

4.81 • 10” 4 

5.03 • 10” 2 

6.73 • 10” 4 


Table 4: Relative errors for the constant coefficient advection diffusion case 
in the /oc-norrri, nonconservative (SL) and conservative (FFSL) scheme, first 
and third order space discretizations. 
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Resolution 

I 2 rel. error (SL) 

I 2 rel. error (FFSL) 

N 

M 

R 

h 

h 

Ro 

R‘2 

50 

50 

0.1 

1.34 • 10" 1 

1.10 • nr 3 

1.30 • 10" 1 

7.31 • 10~ 3 

100 

25 

0.8 

2.10 • 10~ 2 

5.74 • io~ 3 

2.10 • 10" 2 

4.75 • 10~ 3 

100 

100 

0.2 

7.26 • 10~ 2 

4.21 • nr 3 

7.21 • nr 2 

3.52 • lO" 3 

200 

50 

1.6 

1.02 • 10~ 2 

4.68 • io~ 3 

9.57- 10" 3 

3.00 • nr 3 


Tabic 5: Relative errors for the variable coefficient pure diffusion case in the 
/ 2 -norm, nonconservative (SL) and conservative (FFSL) scheme, first and 
third order space discretizations. 

6.2 Linear diffusion with variable coefficients. 

The diffusion equation 

ut = {v(x,t)u x ) x x€[0,L\ 

was then considered on a space interval [0, L\ and time interval [0, T] 
with L = 10 and with T = 4, with time steps defined as At = T/M. 
Periodic boundary conditions were assumed and a Gaussian profile 
centered at L/5 was considered as the initial condition. The diffusivity 
field was given by 


v(x, t) = -1— £(x) sin 

respectively, where £(x) denotes the characteristic function of the in¬ 
terval [0.5L, 0.8L], This choice highlights the possibility to use the 
proposed method seamlessly also with strongly varying diffusion co¬ 
efficients. In this case, no exact solution is available and reference 
solutions were computed using the finite difference method described 
in the previous section with a four times higher spatial resolution, cou¬ 
pled to a high order multistep stiff solver in which a small tolerance 
and maximum time step value were enforced. A plot of the numerical 
solutions obtained at the final time T are displayed in Figures [l]{2j for 
the FFSL and SL method, respectively. 

A more quantitative assessment of the FFSL solution accuracy 
can be gathered from Tables [5} [6j It can be observed that the FFSL 
method is always of slightly higher accuracy than the corresponding 
SL variant at equivalent resolution. In all these computations, FFSL 
maintains mass conservation up to machine accuracy by construction, 
while the average conservation error of the SL method is of the order 
of 10 -3 of the total initial mass. 
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Resolution 

Z 2 rel. error (SL) 

I 2 rel. error (FFSL) 

N 

M 

P 

h 

h 

Ro 

R‘2 

50 

50 

0.1 

1.50 • 10" 1 

i.50 • nr 2 

1.47- 10" 1 

i.28 • nr 2 

100 

25 

0.8 

2.72 • 10" 2 

7.97 • 10- 3 

2.77- 10" 2 

8.21 • nr 3 

100 

100 

0.2 

8.78 • 10" 2 

7.42 • 10" 3 

8.87 • 10 -2 

6.71 • nr 3 

200 

50 

1.6 

1.36 • IQ' 2 

6.95 • 10~ 3 

1.43 • 10~ 2 

5.60 • nr 3 


Tabic 6: Relative errors for the variable coefficient pure diffusion case in the 
Zoo-norm, nonconservative (SL) and conservative (FFSL) scheme, first and 
third order space discretizations. 


6.3 Gas flow in porous media. 

We reconsider the nonlinear example, already treated in [2], of the 
one-dimensional equation of gases in porous media, 

= (23) 


focusing in particular on the so-called Barenblatt-Pattle self-similar 
solutions (see, e.g., hd, which can be written in the form 


u(x,t ) 


(t + t 0 ) k CA 2 


k(m — l)|a;| 2 \ m ~ 1 

2m(t + t 0 ) 2k J , 


(24) 


where to > 0, A is an arbitrary nonzero constant and k = l/{m + 1). 
In order to adapt the FFSL scheme to this case, we recall what has 
been already observed for the nonconservative scheme in [2]: since the 
solution (and hence, the diffusivity) may have a bounded support, a 
straightforward extension of ([b]) in the linearized form 

St = sj2Atu(l[V^}( Xi ±5t)) (25) 


(where in our case v(u) = mu m_1 ) would have multiple solutions for 
Xi out of the support but in its neighbourhood - to ensure a correct 
propagation of the solution, Sf should be defined as the largest of 
such values. The same caution should be applied in extending the 
definition (10). A possible answer in this respect is to compute the 
via (25), then define 5^ = (d~£ + 6^)/2, where, of course, Xk is now a 
cell interface. 

Figure [3] compares the exact and approximate evolution of the 
Barenblatt-Pattle solution for (23), with m = 3, A = 1 and to = 1- 
Approximate solutions have been computed with the nonconservative 
(a) and the conservative (b) scheme at T = 1,4,16 on a mesh com¬ 
posed of 50 nodes with At = 0.05, using respectively cubic interpola¬ 
tion and quadratic reconstruction. Note that the mass conservation 
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Resolution 

I 2 relative error (SL) 

I 2 relative error (FFSL) 

N 

M 

h 

/ 

3 

Ro 

R‘2 

50 

320 

0.316 

8.69- 

10" 2 

0.270 

2.67- 10" 2 

100 

640 

0.212 

4.76- 

I — 1 

O 

1 

to 

0.156 

i.65 • nr 2 

200 

1280 

0.171 

4.84- 

h- 1 

O 

1 

to 

9.85 • 10- 2 

1.09 • 10” 2 

400 

2560 

0.135 

4.49 • 

I — 1 

O 

1 

to 

6.03 • 10~ 2 

2.69 • 10“ 3 

800 

5120 

0.107 

3.22 • 

1— 1 
0 

1 

to 

3.72 • IQ’ 2 

2.79 • 10“ 3 


Table 7: Relative errors for the Barenblatt-Pattle solution in the 2-norm, 
nonconservative (SL) and conservative (FFSL) scheme, first and third order 
space discretizations. 


constraint apparently improves the accuracy of the scheme, especially 
for larger simulation times. This behaviour is confirmed, in terms of 
both absolute accuracy and convergence rate, by Table [7j which shows 
the numerical errors for the two schemes in the 2-norm, at T = 16, 
under a linear refinement law. SL scheme has been tested with linear 
(/i) and cubic ( I 3 ) interpolation, whereas FFSL scheme with piecewise 
constant (Rq) and piecewise quadratic (R 2 ) reconstruction. 


6.4 Variable coefficient case in two space di¬ 
mensions. 

As a two-dimensional test, we consider the equation 


ut = di v(u(x)Vii) 

on 0 = [—3, 3] 2 with periodic boundary conditions. The initial condi¬ 
tion given by the characteristic function of the square £ = [—1.5,1.5] 2 . 
The isotropic diffusivity v is given by 

v{x) = e “ 5 l !C “ x °l 2 , 


where xq = (1.5,—1.5). The diffusion is therefore concentrated in a 
corner on the boundary of the set £. The effect of this diffusion is to 
move mass from the interior of £ to the exterior, in the neighbourhood 
of the point x$. Figure [4] shows the numerical solution along with its 
level curves at T = 2, with a 50 x 50 space grid, Rq reconstruction and 
time step At = 0.05. Note that, despite being conceptually obtained 
by directional splitting, the two-dimensional scheme (21) does not 
suffer from anisotropies induced by grid orientation. 
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7 Conclusions. 


We have extended the SL approach of [2j to achieve a fully conserva¬ 
tive, flux form discretization for linear and nonlinear parabolic prob¬ 
lems. A consistency and stability analysis of the method has also been 
presented, along with a strategy to couple the FFSL discretizations of 
advection and diffusion terms. A number of numerical simulations val¬ 
idate the proposed method, showing that it is equivalent in accuracy 
to its nonconservative variant, while allowing to maintain mass con¬ 
servation at machine accuracy. The proposed method could represent 
an important complement to the many conservative, flux-form semi- 
Lagrangian schemes for advection proposed in the literature. Future 
developments will focus on the application of the proposed approach 
to conservation laws with nonlinear parabolic terms, such as e.g. the 
Richards equation, and to the extension of this technique to high order 
discontinuous finite elements discretizations such as those proposed in 

m, G23- 
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(b) 


Figure 1: One-dimensional case with variable coefficients: numerical solution 
by FFSL method with (a) linear (b) cubic reconstruction. Numerical solution 
is represented by circles, reference solution by continuous line). 
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(b) 


Figure 2: One-dimensional case with variable coefficients: numerical solution 
by SL method with (a) linear (b) cubic reconstruction. Numerical solution 
is represented by circles, reference solution by continuous line). 


22 







Figure 3: Evolution of the exact (continuous line) versus numerical (circles) 
Barenblatt-Pattle solution for T=l,4,16, nonconservative (a) and conserva¬ 
tive (b) scheme. 
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(b) 

Figure 4: Variable isotropic diffusion, graph (a) and level curves (b) of the 
solutions. 
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